/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Implementations of the non-templated ray functions.

// $Id$

#include "ray.h"

/** Scatters a ray in a specified direction. The angles are given in
    the ray frame, in a system where the z-axis is aligned with the
    ray direction and the orthogonal projection of the world z-axis
    and the ray direction is the x-direction. */
void mcrx::ray_base::scatter_in_direction (T_float cos_theta, T_float phi)
{
  assert (cos_theta >= -1);
  assert (cos_theta <= 1);

  const T_float sp = sin(phi);
  const T_float cp = cos(phi);
  const T_float sintheta=sqrt(1-cos_theta*cos_theta);
  const vec3d scatterdir(sintheta*cp,sintheta*sp,cos_theta);

  // Scatterdir is given in a system where the z-axis is aligned with the
  // ray direction. Now we need to rotate this system into the world system
  // We choose the orthogonal projection of the world z-axis as the second
  // axis vector
  const vec3d z(0.,0.,1.);
  vec3d k1(z-direction ()*dot(z,direction ()) );
  // If the ray is nearly aligned with the z-axis this will be a singular
  // system. In that case we just assign scatterdir as the new direction
  const T_float mag_k1 = sqrt (dot (k1, k1)); 
  if(mag_k1 > 1e-5) {
    k1/=mag_k1;
    // And the third axis is the cross product of the two first
    const vec3d k2(cross(direction (),k1));
  
    // Now we perform the matrix multiplication in written-out form  
    const vec3d newdir=
      scatterdir[2]*direction ()+scatterdir[0]*k1+scatterdir[1]*k2;
    
    set_direction (newdir);
  }
  else {
    // Ray parallel to z-axis. This means in this case we use the
    // world y and x axes as the ray x and y axes for the rotation.
    const vec3d newdir=
      scatterdir[2]*direction ()+scatterdir[0]*vec3d(0,1,0)+
      scatterdir[1]*vec3d(1,0,0);
    set_direction (newdir);
  }

  // zero the traversed length and column density
  n_ = 0;
  l_ = 0;
  ++nscat_;
}


mcrx::T_float mcrx::ray_base::intersect_cube(const vec3d& min, const vec3d& max, 
					     bool assume_inside) const
{
  // Calculate lengths until intersections with the cube planes
  const vec3d l1 = (min - position())*inverse_direction();
  const vec3d l2 = (max - position())*inverse_direction();

  // Find out which is entry and exit
  const T_float lx_min = std::min (l1 [0], l2 [0]);
  const T_float lx_max  = std::max (l1 [0], l2 [0]);
  const T_float ly_min = std::min (l1 [1], l2 [1]);
  const T_float ly_max  = std::max (l1 [1], l2 [1]);
  const T_float lz_min = std::min (l1 [2], l2 [2]);
  const T_float lz_max  = std::max (l1 [2], l2 [2]);

  const T_float l_entry = std::max (lx_min, std::max (ly_min, lz_min));
  const T_float l_exit  = std::min (lx_max, std::min (ly_max, lz_max));

  // See if it hits
  if(l_entry > l_exit) // miss
    return blitz::quiet_NaN(T_float());

  if(assume_inside) {
    // if we are on the inside, we return l_exit, but l_entry should
    // not be positive
    assert(l_entry<1e-10);
    return l_exit;
  }
  else {
    // if we are on the outside, we return l_entry, unless both
    // l_entry and l_exit are negative in which case we're heading away
    if(l_entry>=0)
      return l_entry;
    else
      return blitz::quiet_NaN(T_float());
  }
}

